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Abstract 


A  practical  method  for  the  calculation  of  the  spectral 
parameters  of  first-arriving  signals  in  seismology  has  been  the 
object  of  much  theoretical  work  in  recent  years.  The  difficulty 
has  been  in  the  fact  that  the  first  arrival  usually  behaves  as 
an  imperfectly  trapped  mode.  Mathematically,  it  arises  from 
the  contributions  of  branch  line  integrals  and  complex  poles. 
Attempts  to  transform  the  solution  into  a  generalization  of  the 
normal  modes  have  been  a  mathematical  success  only.  Because  of 
the  complexity  of  this  solution,  a  different,  less  elegant 
approach  is  demanded. 

A  practical  technique  is  proposed.  By  a  change  of  variable, 

the  twice  transformed  solution  is  separated  into  a  Droduct  of 

-iuR 

the  form  f(u).e  .  This  can  be  integrated  with  respect  to  the 
phase  variable  u,  using  standard  quadrature  methods,  with  the 
real  part  of  u  changing  most  rapidly  along  the  integration  path. 
By  making  the  frequency  complex,  it  is  possible  to  displace  any 
singularities  away  from  the  vicinity  of  the  contour.  This  gives 
the  spectrum  of  the  signal  as  viewed  through  an  exponentially 
decaying  time  window,  making  it  possible  to  work  with  the  first 
arrival  by  itself. 


Introduction 


The  problem  of  an  impulsive  point  source  or  line  source 
in  a  layered  elastic  medium  is  central  to  the  field  of  seismology 
This  is  due  as  much  to  the  fact  that  many  earth  problems  can  be 
modeled  by  a  series  of  flat  layers  or  spherical  shells,  as  to  the 
ease  with  which  formal  integral  solutions  can  be  written  down. 
Basic  methods  are  discussed  by  Ewing,  Jardetzky,  and  Press  [1957] 
and  algorithms  for  treating  the  complete  n-layered  point  source 
problt  are  presented  by  Harkrider  [1964]. 

These  solutions  are  useful  only  to  the  extent  that  they  can 
be  numerically  evaluated.  For  example,  the  response  to  a  line 
source  (two-dimensional  problem)  has  the  form* 

f(x,z,t)  8  r  [  e  f  cos(kx)  f  dk  du  (1) 

2 71  ;  M«,k) 
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and  the  resoonse  to  a  point  source  (three-dimensional  problem) 
has  the  form+ 


f  (r , z , t ) 
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dk  du 
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The  inner  integrals 

F  (x ,  z , w)  = 


and 


f  cos  (kx)  Si  dk 

J0  6 


(3) 


F(r,z,w)  = 


k  J  (kr)  2.  dk 
>  n  A 


(4) 


give  the  Fourier  transform  of  the  impulse  response.  Since 
observed  seismic  signals  may  be  readily  transformed,  evaluation 


of  (3)  or  (4)  may  be  considered  an  adequate  numerical  realization 
*f  is  any  observable  such  as  a  displacement  or  stress  component 

tn  is  usually  0  or  1,  but  may  be  larger  for  multipole  sources. 

See  Harkrider  [1964]. 
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of  the  theory. 

Evaluation  of  the  integral  has  always  proceeded  by  some  kind 
of  transformation  of  the  problem.  Pekeris  [1948]  deformed  the 
contour  of  integration  in  the  complex  k-plane  in  a  way  which 
separates  the  residue  contributions  due  to  real  poles  of  the 
integrand  (Fig.  la  and  lb) .  All  singularities  lie  on  the  real 
k-axis  if  it  is  reauired  that  all  radicals  of  the  form  v  = 

[k2  -  (w/vn) 2]  1//2  be  taken  with  Re(v)'-  0  on  the  sheet  of  integra¬ 
tion.  Branch  points  at  k  =  w/a  and  k  =  u)/6_  constrain  the 

n  n 

transformed  contour  C'  in  Fig.  lb.  All  poles  lie  to  the  right 
of  the  branch  points  and  to  the  left  of  the  limiting  point  at 
k  =  u)/v  .  ,  where  v_-  is  the  least  velocity  in  any  layer.  An 

u>-k  diagram  (Fig.  2)  summarizes  these  relationships. 

This  transformation  separates  tne  solution  into  a  sum  of 
residue  contributions  and  a  branch  line  integral.  The  former, 
having  phase  velocities  less  than  any  halfspace  velocity,  are  the 
trapped  modes.  Application  of  the  stationary  phase  method  gives 
the  trapped  modes  as  a  function  of  time  in  a  compact  and  physi¬ 
cally  meaningful  form:  the  signal  appears  as  a  superposition  of 
wave  packets  traveling  with  characteristic  group  velocities  of 
the  modes.  This  formalism  is  the  basis  for  all  work  in  surface 
wave  kinematics,  and  has  proved  to  be  a  satisfactory  way  to  bring 
wave  theory  into  seismo logical  practice. 

Description  of  the  early-arriving  signals  from  integrals  (1) 
through  (4)  has  proven  more  difficult.  These  events  consist  of 
the  various  refracted  body  phases  which  travel  at  high  phase 
velocities —  i.  e.  the  contribution  from  the  brancn  line  integral. 
Physically,  the  refraction  arrivals  may  involve  the  same  types  of 
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resonance  which  produce  the  trapped  modes,  but  the  progressive 
loss  of  energy  into  the  halfspace  gives  the  so-called  "leaking 
modes"  instead.  Mathematically,  these  resonances  in  the  refrac¬ 
tion  arrival  come  about  because  complex  poles  lying  on  a  Riemann 
sheet  where  Re(vn)  <0  are  sufficiently  close  to  the  branch  line 
to  affect  the  value  of  the  branch  line  integral  (Fig.  lb) . 

Gilbert  [1964] ,  having  investigated  in  detail  the  behavior  of 
pole  loci  on  the  lower  sheets,  has  suggested  that  his  complex 
ui— k  diagrams  may  be  useful  in  connection  with  frequency /wave- 
number  analysis  of  data. 

Recent  work  has  concentrated  on  expressing  the  branch  line 
integral  in  the  same  formalism  which  applies  to  the  trapped 
modes  (Rosenbaum  [19601,  Phinney  [1961],  Rosenbaum  [1965]). 

This  involves  complicated  transformations  of  the  contour  onto 
the  lower  Riemann  sheets  in  such  a  way  that  the  complex  poles 
are  "picked  up”  as  residue  contributions.  Certain  practical 
problems  are  involved  in  routinely  calculating  these  leaking 
modes.  In  the  first  place,  one  is  still  left  with  line  integrals 
which  must  be  evaluated  numerically;  the  fact  that  they  repre¬ 
sent  only  non-oscillatory  signals  does  not  eliminate  their  im¬ 
portance  in  the  refraction  problem.  More  bothersome  is  the  pro¬ 
blem  of  deciding  which  complex  poles  contribute  to  the  signal; 
relevant  poles  must  satisfy  an  "accessibility"  condition  in  the 
complex  w(k)  plane  before  a  saddle  point  method  may  be  applied. 
These  problems,  as  veil  as  the  sheer  magnitude  of  the  mathema¬ 
tics,  leave  us,  I  think,  with  a  method  which  can  never  be  used 
routinely  for  problems  in  earth  structure  by  other  than  a 
dedicated  specialist. 


(4) 


Direct  Evaluation  of  the  Integrals  on  a  Computer: 

In  this  paper  my  object  is  to  investigate  the  possibility 
of  numerically  evaluating  integrals  such  as  (1)  -  (4)  on  a  digi¬ 
tal  computer.  For  the  sake  of  discussion,  I  will  concentrate  on 
(3) ,  the  Fourier  transform  of  the  line  source  seismogram.  Fortu¬ 
nately,  the  elegance  of  contour  integration  is  still  at  our 
disposal,  making  it  possible  to  consider  several  ays  to  go  about 
the  integration. 

A  minimum  requirement  is  that  the  integrand  must  vary  smooth¬ 
ly  enough  along  the  contour  to  permit  aoproximation  of  the  integral 
by  a  sum.  When  x  is  large,  this  may  be  troublesome,  owing  to 
the  cos{kx)  factor.  A  method  due  to  Filon  [1928]  avoids  this 
problem  by  taking  account  of  the  cosine  function  in  formation  of 
each  of  the  contributions  to  the  sum.  Essentially  a  generali¬ 
zation  of  Simpson’s  rule,  the  Filon  method  requires  only  that  g/A 
vary  slowly  within  the  interval  of  approximation.  By  considering 
only  cases  for  which  both  source  and  receiver  are  near  or  within 
the  layered  part  of  the  structure,  we  have  only  to  worry  about 
the  behavior  of  g/A  in  the  neighborhood  of  the  contour  of  inte¬ 
gration  (the  restriction  on  source  and  receiver  will  be  removed 
later) . 

Evidently  the  contour  pictured  in  l igure  la  is  unsatisfactory 
by  this  standard,  since  it  passes  over  all  the  singularities.  The 
obvious  solution  is  to  transform  to  C1  in  Figure  lb,  and  remove 
the  trapped  mode  poles  from  the  problem,  leaving  a  fairly  well- 
behaved  branch  line  integral.  While  this  may  be  satisfactory 
for  many  values  of  w,  it  leads  to  difficulty  when  the  integral 
must  be  evaluated  for  a  series  of  closely  spaced  frequencies 
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(for  this  is  what  constitutes  a  spectrum) .  In  the  neighborhood 
of  each  cutoff  frequency,  a  normal  mode  pole  passes  through  the 
branch  point  from  the  lov/er  to  the  upper  Riemann  sheet  (Fig.  2). 
Then  g/A  varies  too  rapidly  near  the  branch  point  for  a  numeri¬ 
cal  method  to  be  usable.  Furthermore,  if  this  pole  is  to  the 
.right  of,  but  very  close  to  the  branch  point,  its  effect  v/iil 
have  already  been  subtracted  out  as  a  residue  contribution,  and 
we  will  be  in  the  position  of  adding  in  a  nearly  equal  contribu¬ 
tion  of  opposite  sign.  It  is  not  possible  to  ignore  these  cutoff 
frequencies,  since  they  provide  the  most  important  contributions 
to  the  refraction  arrival.  Further  manipulation  of  C."  would  lead 
to  transformations  on  the  lower  sheet  like  those  which  have  already 
proven  to  be  too  cumbersome. 

We  are  led  back  to  consideration  of  the  situation  in  F;g.  la. 
Two  apparent  ways  to  smooth  the  integrand  along  the  contour  are 
(Fig.  3a  and  3b) :  a)  deform  the  contour  upward  in  the  k-plane 
to  avoid  the  singularities  and  b)  move  the  singularities  into 
the  fourth  quadrant.  It  is  important  to  note  that  the  second  of 
these  methods  changes  the  problem  itself,  since  an  analytic  func¬ 
tion  is  determined  by  its  singularities.  This  is  done  by  making 
«  a  complex  variable  «  =  w-iX;  the  singularities  then  lie  along 
a  straight  line  in  the  fourth  quadrant  which  passes  through  the 
origin.  An  investigation  of  how  (b)  changes  the  problem  will 
show  that  only  the  modified  problem  is  at  all  reasonable  and  why, 
as  a  consequence,  method  (a)  must  fail. 

Consider  the  Fourier  pair  relating  a  signal  f(t)  and  its 
transform  F(w) ,  where  the  bar  over  w  denotes  that  u>  is  considered 
as  a  complex  variable  Z  =  w-iX: 


(6) 


f(t)  =  -L_ 

-0lu)tFU)  dw 

tt>  t  _ 

(5) 

fu)  =  -L- 

/  e“lu)t  f(t)  dt 

(6) 
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If  f(t)  is  identically  zero  prior  to  t=  0,  then  F  is  analytic  in 
the  lower  half  plane,  and  the  path  of  integration  in  (5)  may  be 
run  along  the  line  ImU)  =  -A  =  constant.  The  functional  values 
of  F  along  this  path  are  precisely  the  same  function  implied  if 
we  evaluate  (3)  for  complex  w,  since  we  have  performed  analytic 
continuation  of  functions  identical  on  the  real  axis  A  =  0. 

Calling  this  modified  transform  F^ ,  we  then  get,  by  a  trivial 
change  of  variable: 

[e"xt  f(t)]  =  -L.  I"  elu)t  F.U)  du)  (7) 

/2?  —  A 

F.(u>)  *  JL  /“  e“iu)t[e~At  f  (t)  ]  dt  (8) 

A  /TtT  u°° 

Thus  the  modified  transform  is  a  function  obtained  by  premulti- 
; lying  the  data  f(t)  by  e  Xt  and  Fourier  transforming  in  the  usual 
way.  We  have,  in  effect,  smoothed  the  spectrum.  It  is  signifi¬ 
cant  thrt  this  comes  about  in  an  attempt  to  find  a  practical  way 
to  do  an  integral;  in  retrospect,  it  will  be  apparent  that  an 
unsmoothed  spectrum  will  be  an  exceedingly  ill-behaved  function 
of  frequency  which  is  difficult  to  represent  by  a  table  of  values 
produced  by  a  computer. 

T.ie  impulse  response  of  a  simple  layered  system  persists 
to  times  greatly  in  excess  of  the  slowest  direct  body  wave  arri¬ 
val  time,  due  to  the  partial  trapping  of  energy  bouncing  at  near 
vertical  angles  in  the  waveguide.  A  faithful  representation  of 
the  Fourier  transform  requires  that  a  very  long  time  interval  be 


T 
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analyzed;  in  the  absence  of  smoothing,  this  spectrum  can  vary 
sharply  within  frequency  intervals  equal  to  the  reciprocal  of  the 
time  interval.  Matters  become  most  serious  in  the  rapid  variation 
of  the  phase  with  frequency,  which  leads,  in  the  case  of  a  com¬ 
puter-produced  table,  to  ambiguity  in  the  total  phase  angle. 

Anyone  who  has  computed  "accurate"  Fourier  spectra  from  shorr- 
period  vertical  signals  can  especially  appreciate  the  problems 
involved. 

The  result  is  that,  if  numerical  results  are  an  objective, 
and  one  must  deal  with  discrete  valued  functions  computed  at 
discrete  frequency  spacings,  the  "true"  unsmoothed  spectrum  is 
worthless*.  It  is  too  sensitive  to  small  errors  and  to  details 
of  the  analysis.  In  fact,  the  most  natural  way  to  smooth  is  to 
premultiply  the  signal  by  a  damped  exponential,  since  the  theore¬ 
tical  problem  can  be  solved  by  evaluating  (3)  at  complex  u>.  The 
alternative  method,  shown  in  Fig.  3a,  is  unacceptable,  since  the 
result  must  be  the  unsmoothed  spectrum;  one  would  find  that 


% 
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amplification  of  the  singularities  by  an  exponential  in  the  imag¬ 
inary  part  of  k  would  cancel  out  any  smoothing  of  the  integrand 
obtained  by  moving  away  from  the  singularities. 

It  can  now  be  said  that  we  have  a  wave  theory  for  first 

arrivals.  How  "first"  in  arrival  we  are  talking  about  lies  in 

the  arbitrary  character  of  X.  An  intelligent  decision  requires 

understanding  of  the  problems  of  a  particular  application, 

and  will  depend  on  such  parameters  as  distance,  source  spectrum, 

and  the  amount  of  detail  desired.  In  an  application  with  little 

t  A  possible  exception  would  be  an  expecially  "leaky"  waveguide 
with  a  very  short  response  time,  such  as  a  high  velocity  layer 
on  a  low  velocity  substratum. 
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advance  knowledge  about  the  unknown  earth  structure,  one  wants 
to  use  a  large  A,  to  provide  strong  smoothing  and  a  relatively 
small  number  of  degrees  of  freedom  in  the  data.  When  data  and 
theory  can  be  reconciled  at  this  level,  another  approximation  to 
tne  correct  structure  may  be  studied  by  working  with  a  longer 
time  window  (smaller  A).  While  it  may  be  premature  to  mention 
this,  a  very  desirable  feature  of  this  approach  is  that  the 
inversion  problem  for  body  wave  spectrums  seems  to  have  a  natural 
solution  in  <-erms  of  this  sequence  of  smaller  and  smaller  A. 

A  salient  property  of  the  solution  is  that  it  is  a  single 
function  which  contains  spectral  information  bearing  in  some 
degree  (depending  on  A)  on  all  portions  of  the  signal.  All 
trapped  modes  as  well  as  all  leaking  modes  and  undispersed  body 
phases  are  described;  only  by  a  choice  of  A  can  the  relative 
contribution  of  the  ] ater-arrivals  be  made  small.  This  is  of 
course  a  limitation  in  working  with  the  better  understood  well- 
dispersed  surface  waves,  for  which  group  and  phase  velocity  are 
'rery  useful  concepts.  In  multilayered  systems,  however,  the 
higher  modes  often  have  a  complicated  group  velocity  behavior 
by  wnich  they  add  up  into  the  body  wave  arrivals  from  intermediate 
layers.  Such  a  mode  sum  is  more  conveniently  looked  at  directly, 
through  its  smoothed  spectrum. 
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Further  Practical  Matters: 

We  have  shown  that  numerical  difficulties  associated  with 
evaluation  of  the  line  source  solution  (3)  can  be  avoided  by 
combining  Filon’s  method  with  use  of  a  complex  frequency.  Before 
looking  at  examples,  it  is  necessary  to  cover  some  matters  of 
technique . 

1.  Truncating  the  infinite  integral:  It  is  impractical  to 
evaluate  the  contour  integral  for  large  enough  values  of  k  to 
reduce  the  truncation  error  to  an  acceptable  (1%)  level .  The 
entire  half-plane  to  the  right  of  the  last  singularity  is  analy¬ 
tic,  allowing  (3)  to  be  split  into  posi*'  *e  and  negative 
exponentials : 

,  ,  oo  # 

F  (X,z,u)  =  .5/e  X(g/A)dk  +  . 5/e~lkx (g/A) dk  (9) 

x  x°° 

and  the  respective  contours  to  be  deformed  into  the  upper  and 
lower  half  planes  as  shown  in  figure  4,  with  the  turning  point 
beyond  the  last  singularity.  These  two  line  integrals  converge 
rapidly,  and  can  be  evaluated  by  any  of  the  standard  quadrature 
methods. 

2.  The  high  frequency  problem:  For  a  fixed  choice  of  X  a 
computed  spectrum  becomes  more  and  more  complicated  as  higher 
frequencies  are  considered.  This  is  due  to  the  interaction  of 


many  "modes",  producing  beats  in  the  spectrum.  When  going  to 
higher  frequencies  with  a  fixed  time  window,  one  is  getting 
a  larger  number  of  degrees  of  freedom  associated  with  the  large 
number  of  cycles  in  the  analysis  interval.  In  a  program  aimed 
at  data  analysis  and  interpretation,  there  is  no  reason  why  X 
cannot  be  increased  in  proportion  to  w,  giving  a  time  window 
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which  is  variable,  but  in  fixed  proportion  to  the  period  being 
analyzed.  The  spectrum  will  then  have  the  same  number  of  degrees 
of  freedom  per  frequency  interval  at  all  frequencies;  u>  may  be 
considered  as  varying  along  the  contour  OA  in  figure  5.  The 
constant  X  case  is  contour  BC.  To  differentiate  them  from  the 
raw,  unsmoothed  spectrum  let  BC  define  the  X-spectrum  and  path 
OA  the  ^-spectrum,  since  the  latter  preserves  th,  phase  angle 
of  w.  Subsequent  discussion  about  the  X-spectrum  will  refer  to 
the  ^-spectrum  unless  exception  is  made. 

An  obvious  extension  of  these  notions  is  to  evaluate  the 
spectrum  on  the  contour  OD  or  EF.  These  are,  respectively,  the 
Laplace  transform  of  the  signal  or  the  Laplace  transform  of  the 

iwt 

signal  as  premultiplied  by  e  .  Complete  knowledge  of  either 
of  these  implies  complete  knowledge  of  the  spectrum  on  any  other 
contour,  such  as  OA  or  BC,  since  F  is  analytic  in  the  lower  half 
plane,  and  we  may  go  from  one  contour  to  another  by  analytic 
continuation.  The  difficulty  is  that  analytic  continuation  can¬ 
not  be  performed  with  numerical  data  on  a  contour,  since  we  do 
not  know  all  the  derivatives  of  F  on  that  contour.  This  is 
exactly  the  same  as  attempting  to  solve  an  elliptic  differential 
equation  by  applying  numerical  operators  to  the  data  on  an  ini¬ 
tial  value  line;  the  problem  is  said  to  be  ill-posed.  The 
same  may  be  said  of  attempting  to  perform  an  upward  continuation 
of  data  on  OA  or  BC;  then,  however,  we  have  a  reasonable  amount 
of  information  in  our  data  and  do  not  want  continuation.  In 
contrast,  the  Laplace  spectrum  is  virtually  devoid  of  information 
which  can  be  seen  above  the  numerical  noise.  The  situation  is 
illustrated  in  Fig.  6,  which  show”  a  Laplace  spectrum  and 
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a  X -spectrum. 


3.  Extension  to  the  point  source  problem:  The  essence  of 
Pilots  method  is  that  the  integral  is  divided  into  subintegrals 


of  the  form 

kL+6 

I  cos (kx) (g/A)  dk  ; 

k .  -6 

if  g/A  varies  slowly  enough,  it  can  be  represented  by  a  second 

degree  polynomial  in  the  interval.  Knowing  the  indefinite  inte- 

2 

gral  jcos(kx) [a+bk+ck  ]dk,  we  can  then  write  down  a  quadrature 
formula,  with  6  a  parameter.  The  equivalent  operations  for  the 
point  source  reduce  to  the  problem  of  knowing  the  indefinite 


integral 


ki+6  p 

I (k . , 5 ,n,p)  =  /  J  (kr) (k-k. )  dk 

krs  n  1 


Term  by  term  integration  of  the  power  series  for  the  integrand 
gives  a  power  series  for  the  integral  in  terms  of  5  and  k^. 

By  application  of  the  theory  of  Taylor  series,  this  power  series 
can  be  re-expressed  as  an  expansion  about  some  new  point  k^,  giving 
a  ranidly  converging  series  in  6  and  (k^-kQ) .  As  far  as  numeri- 

j  cal  objectives  are  concerned,  this  integral  may  be  regarded  as 
a  known  function,  which  can  be  generated  on  a  computer  as  readily 
as  the  trigonometric  functions.  For  values  of  kj  larger  than 
about  15,  it  is  desirable  to  have  an  asymptotic  formula  for  the 
integral.  Manipulations  leading  to  such  a  formula  are  straight¬ 
forward  and  will  not  be  discussed  here.  We  can  conclude  that 
the  point  source  solution  requires  a  little  more  programming  than 

;  the  line  source,  but  does  not  differ  in  any  substantive  way. 

Jt 

j* 

4.  Extension  of  the  line  source  problem  to  admit  a  source 
in  the  half space:  This  is  best  handled  by  an  illustrative  case. 

»  If  we  assume  a  line  source  and  receiver  both  located  in  the 


(1’) 


half space  (figure  7),  then  the  response  has  the  form 


1 

—  co 


f  -i(kx 

*  A 


+  v  Z) 
a2 


q(w,k)  dk 


(12) 


where  =  .IT2  -  k ^  ;  K  =  u/a_  '  Z 

ex  2 

We  want,  to  put  this  in  the  form 


=  z  +  h 


-  iuR 


F  =  /e  q  (oj,k  (u) )  du 


(13) 


and  integrate  along  a  path  such  that  either  Im(u)  =  constant  or 
Re(u)  a  constant.  The  methods  described  earlier  in  this  paper 
will  then  be  applicable. 

Set  x  =  R  cos  0  k  =  K  cos  n 

I 

Z  =  R  sin  8  va2  =  K  sin  n 


Then  kx  +  va2Z  =  KR  cos(n  -  0) 

I 

or  u  =  K  cos(n  -  0)  -  k  cos  0  +  ua2sin  0  (14) 

We  can  solve  for  k(u)  in  this  special  case: 

k  =  u  cos  0  +  /K7  -—u7’  sin  0  (15) 


with 


dk  = 
du 


cos  0 


u 


I  sin  e 


(16) 


We  are  in  a  position  to  make  use  of  (13) .  The  mapping  of  the 
real  k-axis  onto  the  u-plane  is  described  by  (14)  and  shown  in 
figures  8a  and  8b.  Point  A,  at  u  =  K  sin  0,  is  the  map  of  k  =  0; 
for  small  k,  u  is  real  and  increasing,  as  can  be  seen  from  the 
derivative 


du  _ 
dk 


cos  0 


k 

/K^-k7 


sin 

' 


0 


U 


(13) 


which  va?  xshes  at  k  =  K  cos  9,  the  well  known  stationary  phase 
point,  wher3  8  =  n  .  Point  S,  at  u  =  K,  is  the  map  of  this 
stationary  phase  point,  and  represents  a  singularity  where  the 
mapped  contour  reverses  direction.  This  singularity  (an  addition¬ 
al  branch  point  in  the  u-plane)  will  provide  an  important  part 
of  the  value  of  the  integral.  Point  B,  at  k=K,  u=K  cos  0, 
is  the  singularity  where  du/dk  is  undefined,  and  the  contour 
moves  into  the  complex  u-plane.  The  ambiguity  in  the  square 

root  is  resolved  by  the  original  condition  that  Re(va2)  > 

•  » 

which,  in  light  of  va2  =  iva2,  requires  that  Im(va2)  <  0.  The 
remainder  of  the  contour  is  a  limb  of  a  hyperbola  in  the  fourth 
quadrant,  asymptotic  to  the  line  -Im(u)/Re(u)  =  tan  8.  The 
remaining  singularities  of  the  integrand  are  folded  into  the 
interval  ASB  on  the  real  axis. 

For  reasons  discussed  earlier,  the  frequency  must  be  taken 
with  a  negative  imaginary  part.  The  mapping  of  the  real  k-axis 
is  shown  in  figure  8b.  The  points  A,  B,  and  S  have  moved  into  the 
fourth  quadrant,  but  are  defined  as  in  the  previous  paragraph. 

All  other  singularities  lie  on  the  line  segment  ASB.  The  contour 
AC  is  the  image  of  the  real  k-axis.  AC  can  be  deformed  into  the 
path  AHJ  without  crossing  any  singularities;  the  integration  from 
A  to  H  is  done  by  Filon's  method  and  the  rapidly  converging  inte¬ 
gral  from  H  to  J  by  any  common  quadrature  formula.  This  final 
transformation  demonstrates  an  important  physical  result —  that 
the  predominant  contributions  to  the  integral  come  from  the 
singularities  between  A  and  S,  including  the  stationary  phase 
contribution  itself.  Singularities  lying  on  the  (folded)  segment 
SB  will  not  contribute  significantly.  In  the  original  variables, 


(14) 


this  means  that  wave  numbers  greater  than  the  stationary  phase 
value  do  not  contribute,  or  that  phase  velocities  less  than  the 
critical  phase  velocity  do  not  contribute. 

To  carry  out  the  integration  it  is  necessary  to  know  the 
complex  values  of  k  corresponding  to  points  on  the  contour  AHJ . 
For  the  problem  just  described,  a  solution  (15)  can  be  written 
down.  A  more  general  case  is  the  partial  ray  integral 


F*  - 


-i(kx  + 


V  h  + 

a  2 


V  6  2  2  J 


q(u,k)  dk 


(17) 


which  describes  an  upgoing  P-wave  coupled  to  the  layer,  coupled 
in  turn  to  a  downgoing  S-wave.  The  defining  equation  for  the 
phase  function  u  becomes 


uR 


kx  +  h  + 

a2 


(18) 


where  R  is  some  characteristic  distance  which  we  would  normally 
take  as  the  optical  distance  for  the  PS  reflection.  Inversion 
to  give  k  as  a  function  of  u  must  be  done  numerically,  since  (18) 
cannot  be  solved  algebraically.  With  a  slight  change  in  notation, 
(18)  is  the  same  equation  which  must  be  solved  to  get  the  Cagniard- 
deHoop  path  used  in  the  Cagniard  method.  None  of  the  further 
elegance  of  this  technique  is  of  use  in  boundary  value  problems 
involving  more  than  one  interface. 

5.  Combining  the  previous  two  problems —  a  poirt  source  in 


the  half space:  The  method  must  embody  features  from  the  preceding 
two  sections;  no  attempt  is  made  in  this  paper  to  do  the  problem. 


(15) 


Later  arrivals,  filtering,  and  the  refraction  problem: 

The  numerical  A-spectrum  described  here  most  naturally  is 
applied  to  the  first  arrival —  invariably  some  sort  of  P-wave. 

It  is  of  practical  interest  to  study  the  extent  to  which  one  can 
look  at  later-arriving  body  waves. 

1.  Partial  ray  expansions:  Mu>,k)  contains  sums  of  expon¬ 
entials  in  the  various  layer  parameters.  Various  binomial  expan¬ 
sions  of  g/A  may  be  formed,  each  corresponding  to  a  particular 
way  of  partially  decomposing  the  solution  into  rays.  An  example 
of  one  term  in  such  an  expansion  is  (17),  where  we  can  identify 
exponential  phase  factors  describing  ray  propagation  and  a  gener¬ 
alized  reflection  coefficient  giving  the  complete  response  of 
the  layers  involved.  Many  of  the  later-arriving  body  phases  can 
be  established  as  first-arrivals  of  a  partial  ray  integral.  One 
must  be  reasonably  sure,  when  looking  at  data,  however,  that  the 
body  phase  being  studied  is  not  contaminated  by  earlier-arriving 
energy.  Somewhat  more  limiting  is  the  fact  that  each  partial  ray 
gives  rise  to  more  than  one  refraction  arrival  in  addition  to  the 
principal  optical  ray  event.  These  refractions  cannot  be  separa¬ 
ted  by  further  decomposition  of  the  integrand.  The  situation  is 
illustrated  in  figure  9.  The  partial  ray  is  defined  by  the  geo¬ 
metrical  path  PS  (and  has  the  same  form  a;;  in  (17)).  Evaluation 
of  the  integral  and  formation  of  the  impulse  response  would  give 
a  seismogram  containing: 

a)  the  first  ‘ival,  a  refracted  event  PpS 

b)  a  refracted  event  PsS 

c)  the  reflection  PS 

No  manipulation  of  partial  rays  can  separate  these  three  events. 


(16) 

In  addition,  the  refraction  PpP  arises  from  a  different  partial 
ray  entirely,  but  may  arrive  so  close  to  PpS  that  separation  of 
the  phases  is  not  possible.  It  would  be  necessary  in  this  case 
to  compute  at  least  these  two  partial  rays. 

2.  Filtering  in  the  w-k  plane:  From  our  understanding  of 
the  kinematics  of  propagation  of  body  phases,  we  may  discriminate 
against  the  first  arrival  on  the  basis  of  its  phase  velocity,  in 
effect  transforming  a  later  arrival  into  the  first  arrival.  Since 
filtering  can  be  done  by  simple  multiplication  in  the  variables  u> 
and  k,  multiplication  of  the  integrand  by  a  velocity  filter  can 

be  used  to  reject  higher  phase  velocities  associated  with  first 
arrivals.  Normal  care  in  the  design  of  the  filter  is  still  re¬ 
quired;  a  conservative  procedure  would  be  to  use  only  filters 
with  physical  analogs  in  simple  array  processing. 

3.  Windowing  vy  convolution:  It  has  been  shown  that  use 

of  the  complex  frequency  variable  provides  the  effect  of  an  expo¬ 
nential  window  in  the  time  domain.  It  would  be  desirable  to  have 
more  general  windows  available,  such  as  a  lagged  exponential, 
which  could  be  used  t  elect  later-arriving  body  waves.  I  can 
see  no  other  way  of  doing  this  except  directly —  by  applying 
a  convolution  operation  in  the  frequency  domain.  This  adds 
another  integration,  which  would  be  as  much  work  as  computing 
the  impulse  response  directly  by  Fourier  synthesis. 

4.  Crossover  of  refraction  lines  from  different  layers: 

An  obvious  application  of  the  A-spectrum  is  to  first  arrivals  at 
intermediate  refraction  distances  where  two  or  more  refraction 
lines  have  nearly  identical  arrival  times.  When  a  refraction, 
such  as  line  2  in  figure  9,  does  not  appear  as  a  first  arrival, 
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its  identification  from  later  arrivals  at  large  distances  can  be 
a  serious  problem.  Even  if  it  appears  as  a  first  arrival  in  a 
short  distance  interval,  the  distribution  of  geophones  may  be 
such  that  it  can  only  be  seen  at  one  station,  in  which  case  the 
slope,  or  even  the  identification  of  the  line  is  in  doubt. 

Since  the  A-spectrum  contains  all  the  information  inherent 
in  the  wave  solution  for  the  time  window  and  frequency  band  under 
consideration,  it  is  appropriate  to  work  with  it  at  several  dis¬ 
tances  near  the  crossover.  The  result  is  an  extended  type  of 
refraction  plot:  a  A-spectrum  versus  distance,  which  constitutes 
a  surface  or  a  nest  of  curves.  It  would  be  the  object  of  the 
refraction  study  to  adjust  the  model  parameters  to  bring  the 
theoretical  and  experimental  spectrum-distance  data  into  agreement. 

This  is  one  of  a  class  of  problems  in  which  geometrical 
optics  gives  more  than  one  arrival  line  with  nearly  the  same 
travel  time.  The  neighborhood  of  a  cusp  in  a  travel  time  curve 
is  of  this  nature.  The  multiple  branches  of  the  PKP  phase  which 

i 

grazes  the  inner  core  are  an  example  in  a  spherical  geometry. 

I 

To  interpret  the  data  in  terms  of  ray  theory,  one  must  be  able 

! 

| 

to  distinguish  separate  arrivals,  no  matter  now  close  they  are  in 
time,  and  to  identify  the  point  of  a  cusp.  Neither  is  possible 
with  seismic  data  recorded  on  an  instrument  with  a  finite  pass- 
band.  In  practice,  conventional  short  period  seismometers  come 
close  enough  in  the  PKP  example  to  show  the  existence  of  several 
branches;  it  cannot  be  shown,  however,  whether  other  branches 
have  been  missed,  or  whether  the  identification  of  the  cusps  is 


correct.  This  latter  problem  comes  about  because  of  the  existence 
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The  principle  of  a  X-spectrum-distance  profile  can  be  used  to 
resolve  problems  of  this  nature.  It  will  be  possible  to  use  con¬ 
ventional  long-period  seismometers  for  wave  theory  interpretations 
of  earth  structure  at  certain  depths  where  diffraction  effects 
are  important:  the  base  of  the  upper  mantle  at  4^0+  km;  the  core¬ 
mantle  boundary;  the  inner  core  boundary.  This  paper  does  not 
describe  any  calculations  of  che  spectrum  in  a  spherical  geometry. 
There  are,  however,  no  real  obstacles  to  doing  so. 


I 
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Calculation  of  phase  velocity,  attenuation,  and  group  velocity; 

Phase  velocity  and  attenuation  are  basic  notions  in  the  kine¬ 
matics  of  traveling  waves,  and  are  easily  applied  in  a  description 
of  trapped  modes.  When  studying  arbitrary  transients  containing 
various  types  of  waves,  these  variables  must  be  regarded  merely 
as  derived  quantities,  defined  by  certain  operations  with  the  X- 
spectrum  of  the  signal. 


Writina  the  soectrum  F ,  =  A  e‘ 


(19) 


we  define  the  phase  velocity  c,  by 

0) 

c  'if!1 

and  the  attenuation  coefficient^  by 


(20) 


1  3A 
Y  ■  ■  K  7x 


(21) 


By  differentiation  of  (19) : 


♦  Aifte1* 


(22) 


and  division  by  (19) : 


1  3F. 

p  3x 
P  X 


_  1_  3^  ^  3  $ 

A  3x  Tx 


(23) 


we  get: 


c  =  _  1  3F 

Im(Fx  3x} 


(24) 


and: 


Y  =  “Re (i  |£) 
Fx3x 


(25) 


d  r 

form  -^X  from  the  theoretical  solution,  differentiate  (3): 


|Ix  = 

3x 


-  f  k  sin(kx) (g/A)  dk 

J  r\ 


(20) 


From  (26)  we  can  compu_e  by  Filon's  method  at  the  same  time 

hat  we  compute  F^.  The  theoretical  phase  velocity  end  attenua- 
cion  follow  from  (24)  and  (25) .  It  will  become  apparent  from 
the  numerical  examples  to  follow  that  these  functions  are  not 
always  very  diagnostic  when  several  modes  or  rays  contribute  to 


the  A-soectrum. 


The  taking  of  a  derivative  by  means  of  an  analytic  relation 
such  as  (26)  does  not  eliminate  the  amplification  of  error  which 
occurs  in  the  derivative  operation.  It  is  found  that,  in  general, 
calculations  of  c  and  y  have  about  one  less  decimal  place  of 
accuracy  than  do  the  calculations  of  F  itself. 

A 

A  derived  group  velocity,  U,  can  be  defined  by 


U  =  J 

t  i! 


(27) 


This  is  based  on  the  fact  that  the  impulse  response  integral 


f  . ,  .  i  (wt  +  $)  „ 
f(t)  =  J  A  (id)  e  du) 


(28) 


receives  its  greates  contribution  at  time  t  from  the  frequency 
for  which  -jj—Cwt  +  $)  =  0*  Equation  (27)  follows.  Like  phase 
velocity,  U  is  of  limited  usefulness  in  looking  at  an  arbitrary 
transient.  Further  difficulty  lies  in  the  necessity  of  either 
programming  the  analytic  form  of  |-^  or  obtaining  the  derivative 
by  simple  differencing.  An  illustrative  example  will  use  the 


latter  method. 
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Numerical  examples: 

Calculations  are  shown  for  the  line  source  response  for 
the  liquid  layer/liquid  halfspace  and  liquid  layer/solid  halfspace 
problems.  Parameters  are  listed  in  Table  1. 


.Model 

a 

6 

P 

1 

2 

? 

2 

LPLTE1 

_ 

_ 

LGL002 

3.0 

- 

2.5 

LOS002 

3.0 

1.732 

2.5 

LOS010 

3.0 

.95 

2.5 

LOS011 

3.0 

.40 

2.5 

LOS 012 

3.0 

.70 

2.5 

LOS013 

3.0 

1.30 

2.5 

Table  1.  Layered  models  used  for  calculations.  For  all  models 
_ =  1.0,  a.  =  1.0,  pt  -  1.0,  for  normalization . _ 

Expressions  for  g/A  are  in  Ewing,  Jardetzky,  and  Press  [1957] , 
chapter  4. 

Model  LPLTE1  is  unique  in  having  a  residue  series  solution 

| 

which  is  complete,  and  which  is  unencumbered  by  a  branch  line 
integral.  To  provide  a  check  of  the  numerical  integration 

program,  the  spectrum  for  this  model  was  obtained  in  two  different 

1  # 

ways:  first  by  substituting  p 2  =  0  in  the  liquid/liquid  kernel 

I 

g/A,  and  then  1/  using  the  residue  series. 

If  an  acceptable  numerical  noise  level  is  1  or  2  in  the  third 
decimal,  then  it  is  found  that  Ak  =  .04  is  a  small  enough  spacing 
provided  that  A  ^  .15  and  x  4  20.  On  the  IBM  7094,  calculation 
of  100  spectral  values  takes  about  5  minutes  for  the  liquid/liquid 
problem  and  about  7  minutes  for  the  liquid/solid  problem.  When 
computing  a  A-spectrum  for  a  set  of  different  x,  however,  the 

|| 

same  values  of  g/A  are  generated  in  each  pass.  By  saving  these 


|  pBHiii!Hpffl!lllM!l»|*B'  . . 


(22) 


values  in  a  high  speed  backup  storage  (disk  file)  and  reusing 
them  in  the  second  and  succeeding  passes,  the  time  requirement  is 
reduced  to  about  40  seconds  per  case  for  all  but  the  first. 

While  the  limitation  to  a  >  .15  is  not  serious,  the  restric¬ 
tion  that  x  not  be  too  large  is  annoying.  This  comes  about  from 
accumulation  of  roundoff  error  in  forming  the  sum  of  contributions 
frcm  the  various  k-intervals.  As  x  increases,  one  is  adding  up 
the  same  table  of  g/A  values  with  different  weights;  for  larger 
x,  the  values  tend  to  cancel,  and  roundoff  becomes  a  problem. 

To  a  limited  extent,  this  can  be  overcome  by  grouping  terms  or 
by  selecting  xAk  to  be  a  multiple  of  2 n  aid  grouping.  A  more 
satisfactory  resolution  of  the  problem  is  deferred  at  this  time. 

The  notion  of  a  profile  of  A-spectrum  versus  distance  is 
illustrated  in  figure  11,  which  shows  model  LDL002 ,  with  A  =  .15, 
h  =  1.3,  and  z  =  0.1.  For  a  liquid/liquid  case,  the  early  arriv¬ 
ing  signal  is  very  nearly  described  as  a  sun  of  dispersed  normal 
modes  (Pekeris  [1948]),  with  the  branch  line  integral  contribu¬ 
ting  a  relatively  weak  signal.  Our  curves  should  reflect,  there¬ 
fore,  the  mode  character  of  the  response.  Most  of  the  energy  is 
concentrated  around  the  theoretical  cutoff  frequencies  of  the 
first  two  modes  at  .26  and  .78,  due  to  the  source  being  in  the 
half space.  The  irregular  behavior  (for  x  =  7  and  10)  at  the 
higher  frequencies  is  due  to  the  effect  of  mode  mixing. 

Figure  12  gives  phase  velocities  and  a  single  instance  of 
group  velocity  for  this  model.  The  mode-like  character  of  the 
signal  is  demonstrated  by  the  invariance  of  c  with  x.  The  effect 
of  the  branch  line  integral  is  automatically  included  at  frequen¬ 
cies  near  cutoff,  causing  c  and  U  to  behave  differently  from  c 
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given  by  mode  theory.  The  variability  of  c  in  the  vicinity  of  the 
second  node  cutoff  (due  to  node  interference)  demonstrates  the 
kind  of  limited  service  that  can  be  expected  of  this  parameter  at 
higher  frequencies. 

A  comparable  example  for  the  liquid/solid  problem  LOS002 
appears  in  figure  13.  The  irregular  structure  of  the  amplitude 
function  is  a  result  of  mixing  the  first-arriving  PL  mode  with 
the  trapped  modes.  Even  though  A  =  . 3  in  this  case,  giving  a 
shorter  time  window,  the  trapped  nodes  still  arrive  in  the  time 
covered  by  the  window  except  at  the  greatest  distances. 

The  effect  of  the  basement  shear  velocity  at  a  given  distance 
is  shown  by  figure  14.  In  all  other  respects,  this  analysis  has 
the  same  parameters  as  the  liquid/liquid  case  in  figure  12.  By 
letting  B2  become  small  we  have  a  check  on  the  liquid/solid 
calculation.  When  it  is  as  small  as  0.4,  the  agreement  with  the 
liquid/liquid  result  is  good;  for  values  greater  than  1.0,  the 
spectrum  is  mainly  determined  by  the  existence  of  the  trapped 
shear  modes.  At  the  lower  shear  velocities,  the  shear  and  surface 
wave  energy  arrives  so  late  that  it  is  excluded  by  the  exponential 
time  window,  and  we  see  only  the  effect  of  the  P  modes.  The  phase 
velocity  (figure  15) ,  being  a  derivative  quantity,  is  rather  ill- 
behaved  by  comparison  with  the  spectrum  itself,  and  the  agreement 
with  the  liqui J/liquid  case  is  achieved  only  at  the  lowest  values 
of  6 2 .  Discontinuities  in  c  appear  near  minima  in  the  amplitude 
spectrum. 

My  viewpoint  in  this  paper  has  been  that  the  A-spectrum  is 
an  appropriate  function  to  work  with  in  data  analysis.  The  pre¬ 
ceding  examples  demonstrate  that  its  behavior  may  be  sufficiently 
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strange  to  most  workers  that  a  demonstration  is  in  order.  The 
A-spectrum  has  been  computed  for  LOSOC2  at  x  =  20.  for  an  increas¬ 
ing  sequence  of  A  (figure  16).  Fourier  synthesis  has  been  applied 
to  each  spectrum  to  yield  a  theoretical  seismogram--  as  modified 
by  the  exponential  window.  To  reduce  the  effect  of  chopping  the 
A-spectrum  at  a  particular  frequency,  a  realizable  low-pass 
filter  has  been  applied  before  synthesis.  The  results  appear 
in  figure  17.  The  creation  of  a  theoretical  seismogram  appears 
to  be  the  only  process  in  which  the  A-spectrum  is  preferable  to 
the  ({/-spectrum  (which  has  the  virtue  of  smoothness  at  high  fre¬ 
quency)  .  If  Fourier  synthesis  is  based  on  a  ({/-spectrum,  the 
theoretical  time  function  is  a  peculiar  modification  of  the 
full  seismogram,  in  which  each  frequency  has  been  premultiplied 
by  a  different  time  window.  This  time  function  would  appear  as 
though  it  had  been  generated  by  a  set  of  ccnstant-Q  delay  lines, 
in  which  the  attenuation  coefficient  is  proportional  to  frequency. 
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Conclusion : 

It  has  been  shown  that  the  requirement  of  nume“ically  realiz¬ 
ing  the  theoretical  response  of  a  layered  medium  leads  naturally 
to  the  conclusion  that  signals  must  be  viewed  through  a  decaying 
exponential  time  window.  The  ultimate  object  of  such  work  is 
to  be  able  to  use  comparisons  between  theory  and  data  to  deter¬ 
mine  velocities  and  densities  in  the  earth.  At  the  present  time, 
however,  data  processing  technique  is  considerably  more  advanced 
than  the  theoretical  work,  and  a  considerable  a'  ount  of  computer 
programming  must  be  done  to  close  the  gap.  This  includes  not 
only  complete  realizations  of  the  methods  described  in  the  text 
(or  alternate  methods) ,  but  a  realization  of  available  formal 
schemes  for  obtaining  the  response  function  in  multilayered  media. 
From  a  practical  point  of  view,  the  Thomson-Haskell  matrix  method 
should  be  supplemented  or  replaced  by  methods  based  on  a  continu¬ 
ous  velocity-depth  function  which  will  be  faster  by  virtue  of 
their  dependence  on  fewer  independent  parameters.  For  example, 
the  base  computing  time  for  one  spectrum  is  5  minutes  in  the  most 
elementary  case;  if  20  solid  layers  were  required  to  model  an 
earth  structure,  the  time  could  go  over  200  minutes  per  curve — 
for  the  same  amount  of  output  information.  A  compensating 
factor  will  be  the  availability  in  the  near  future  of  computers 
which  can  provide  10-100  times  the  computing  speed  of  present 
models . 

Despite  this  present  lack  of  theoretical  capability  to  match 


i  with  data  analyses,  it  is  clear  that  work  in  data  processing  of 
body  waves  must  be  affected  by  the  importance  of  the  exponential 

Hi 

window.  Demonstrations  must  be  made  of  the  way  in  which  the 
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A-spectrum  varies  with  distance  and  choice  of  smoothing  parameter. 
In  cases  where  the  source  is  controlled  or  sufficiently  kn  ;wn , 
it  is  necessary  to  establish  the  repeatibility  of  the  result  over 
a  given  path,  and  to  determine  the  lelation  between  A  and  the 
repeatibility.  In  the  area  of  data  analysis  techniques,  it  does 
not  appear  that  work  with  sets  of  analysis  functions  more  sophis¬ 
ticated  than  the  damped  exponentials  can  be  meaningfully  related 
to  wave  propagation  theory  .  This  comes  back  to  the  elementary 
fact  that  the  complex  exponentials  are  the  correct  choice  of 
functions  for  separation  of  the  wave  equation. 
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Appendix 


1.  Conventional  geometry  of  a  layered  system: 
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2.  Notation: 

A  amplitude  of  the  spectrum 

c  phase  velocity 

f  some  field  variable  in  the  time  domain 

F  some  field  variable  in  the  frequency  domain 

h  depth  of  source 

thickness  of  ith  layer 
k  horizontal  wave  number 

n  used  as  a  subscript  referring  to  the  hwlfspace 

r  radial  coordinate  of  receiver  in  cylindrical  coordinates 

t  time 

u  phase  variable  (equation  14) 

U  group  velocity 

v  denotes  a  body  velocity,  either  a  or  8 

x  horizontal  distance  from  source  to  receiver  (line  source 

problem) 


2.  Notation  (continued): 
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6i 
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a1 


6i 
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compressional  velocity  in  ith  layer 
shear  velocity  in  ith  layer 
attenuation  coefficient 
ray  angle  (equations  13-14) 
wave  angle  (equations  13-14) 

(negative  )  imaginary  part  of  complex  frequency 

_  I 

/ky-  ( w/a  )  *  =  iv 

l  ai 

/k^- (u/6 . ) 7  =  iv 

i  6i 

density  in  ith  layer 

phase  angle  of  spectrum 

minus  phase  angle  of  u 

angular  frequency 

complex  frequency  variable  =  w  -  iX 


3.  Fi Ion's  method: 

This  is  a  method  of  numerically  integrating  trigonometric 

integrals  of  the  form 

b 

T  =  /  cos(xp)  f (p)  dp 
a 

where  x  may  be  large.  Let  the  range  be  divided  into  2n  equal 
parts  with  an  interval  h  such  that 
b  =  a  +  2nh 

After  representing  f(p)  over  the  interval  [a+(s-l)h  ,  a+(s+l)h] 
by  a  second  degree  polynomial  whose  coefficients  can  be  determined, 
we  get,  after  some  manipulation 

I  =  h(a{f (b) sin(xb)  -  f(a)sin(xa)}  +  8C2S  +  yC2S-]J 
where,  if  0  =  xh, 

a  =  6~  3 { 6 2  +  0sin9  cos0  -  2sin  28) 


8  =  6  3 { 6  { 1  +  cos2 6]  -  2  sin6  cos9} 
y  =  4  0'3{  sinO  -  0  cos©} 

and 

n  i 

C 2 g  =  1  f (a+2 jh) . cos [x (a+2 jh) ]  -  ^[f (b) cos (bh) +f (a) cos (ah) ] 
j=0  2 

n 

C9  .  =  l  f [a+(2j-l)n] .cos{x[a+(2j-l)h] } 
j*l 

A  complete  description  of  the  method  is  given  by  Tranter  [1956] , 
following  the  original  paper  by  Filon  [1928] . 


Figure  Captions 


la,  k-plane  for  integration  of  (3),  and  initial  contour, 

showing  singularities  on  the  real  axis. 

lb.  k-plane  for  integration  of  (3) ,  with  the  contour  deformed 

into  path  C'  around  the  branch  cut. 

2.  uj-k  diagram  for  singularities  of  the  integrand  on  the  top 
sheet . 

3a.  Numerical  integration  by  deforming  the  contour  away  from 
the  singularities. 

3b.  Numerical  integration  by  moving  the  singularities  into 
the  fourth  quadrant. 

4.  Truncation  of  the  integral  by  splitting  into  exponentials 

and  deforming  the  respective  contours  into  the  upper 
and  lower  half-planes  respectively. 

5.  w-plane.  Examples  of  lines  along  which  the  spectrum  may  be 

computed . 

6.  Comparison  of  computed  Laplace  and  A -spectrum  for  model  LOS002. 

Curves  are  separately  normalized. 

7.  Geometry  for  a  source  and  receiver  in  the  half space. 

8a.  u-plane,  showing  map  of  the  real  k-axis,  using  (14),  when 
a)  is  real. 

8b.  u-plane,  showing  map  of  the  real  k-axis,  using  (14)  ,  when 
Im(io)  is  negative.  The  deformed  contour  AHJ  is  used 
in  the  numerical  evaluation  of  (13) . 

9.  Geometry  for  a  PS  partial  ray  integral  in  a  three-layer  model. 

10.  Illustrative  travel-time  curve  for  a  three-layer  model  with 
a  masked  refraction  arrival. 

11.  Model  LOL002:  A-spectrum  as  a  function  of  distance  for 
h=  1.3,  z=  0.1,  A=  .15,  due  to  an  impulsive  line  pressure 
source . 

12.  Model  LOL002:  Phase  velocity  spectrum  as  a  function  of  distance 
and  a  single  example  of  group  velocity,  obtained  from  the 
A-spectrum  of  Figure  11. 

13.  Model  LOS002:  A-spectrum  as  a  function  of  distance  for 
h=  0.9,  z=  0.1,  A  =  .30,  due  to  an  impulsive  line  pressure 
source. 


Figure  captions  (continued) 


14.  Depend  , nee  of  X-spectrum  on  e2»  f°r  a  sequence  of  models 
with  decreasing  shear  velocity  in  the  half space,  at 

x=  10.,  h=  1.3,  z=  0.1,  X  =  .15.  The  limiting  case  of  B2=  0 
is  computed  using  a  liquid/liquid  program. 

15.  Phase  velocity  spectrum  for  the  models  computed  in  Figure  14. 

16.  Model  LOS002:  X-spectrum  as  a  function  of  X  for  x=  20., 
h=  1.3,  z=  0.1  due  to  an  impulsive  line  pressure  source. 

17.  Model  LOSO"  :  Impulse  response  obtained  by  Fourier  synthesis 
of  spectrum,  in  Figure  16. 
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Figures  8a  and  8b 
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